classdef BPTES < handle
    % 理想气体，Re数变化较小
    % 液态储热介质，不妨就采用天野2023论文中介质
    % 待优化的自由度：T_Charge、T_Discharge以及充放电的压比（认为相同）
    properties
        timeStep = 3600 % 1小时
        tol = 1e-2
        maxIter = 1e4

        powerNominalCharge % 额定功率，充电功率
        powerCharge
        powerDischarge

        capacityBound = [0.1, 0.9] % 热罐能量上下
        capacityEnergy % 此应为T_hot即热罐中能量0~100%之差
        capacityMassGain % 控制热罐与冷罐容量不等的变量
        T1SaltHot = 500 + 273.15 % 单位K
        T2SaltHot = 300 + 273.15
        T1SaltCold = 0 + 273.15
        T2SaltCold = -80 + 273.15
        % 熔融盐的比热容 [J/kg-K]，=1443+0.172*℃，此处记录两个区段均值
        c_pSaltHot = 1511.8
        c_pSaltCold = 1530 % 实际为Hexane的，但此处统称为盐

        % 动力循环工质参数
        c_pCycle = 1060 % 动力循环工质比热容，此处采用空气J/kg·K
        gamma = 1.4 % 动力循环工质比热比
        pressureRatio % 压比

        T1Charge % 待优化自由度，受pressureMin、热罐温度限制，单位K
        % T2Charge
        % T3Charge % 与回热有关 3→2
        T4Charge % 待优化自由度，受冷罐温度限制
        % T5Charge
        % T6Charge % 与回热有关 6→5

        T1Discharge % 待优化自由度，受pressureMin、热罐温度限制
        % T2Discharge
        % T3Discharge % 回热 与T6对应，取决于换热器
        T4Discharge % 待优化自由度，受冷罐温度限制
        % T5Discharge
        % T6Discharge % 回热 T5升温到T6，取决于换热器


        % 设备效率
        efficiencyCompressor = 0.87 % 等熵效率
        efficiencyTurbine = 0.92 
        efficiencyEl = 0.95
        efficiencyGearBox = 0.97

        pinchpointHot = 2 % 2K  
        % T1Charge-T1SaltHot≥pinchpointHot  T6Charge-T2SaltHot≥pinchpointHot
        pinchpointCold = 2 % 2K 
        % T1SaltCold-T3charge≥pinchpointCold  T2SaltCold-T4charge≥pinchpointCold
        pinchpointReg = 2 % 2K
        % T6Charge-T2Charge≥pinchpointReg  T5Charge-T3Charge≥pinchpointReg
        temCompressorMax = 580+ 273.15 % 580℃
        % T1Charge≤temCompressorMax
        pressureMin = 2.5 
        pressureAtm = 1.01325 % 大气压，单位MPa

        % 动态参数
        massSaltHot1  % 热罐1质量
        massSaltHot2 
        massSaltCold1 
        massSaltCold2 
        % massRateCycCharge % 由额定功率推导
        % massRateCycDischarge % 由关联式得到

        % 换热器
        % 换热器Q = f*UA*LMTD，其中LMTD为对数平均温差，f为相对于逆流的修正系数（空气冷却=0.8，其余=0.95）
        % 以放电过程为例，UA为自由度，而需要计算的只有T6，由对称得到T3，由UA与m流率得到T6.这里不用管回热度，因为计算T6本质上就是算回热度

        f_Reg = 0.95
        UA_Reg
        UA_Hot
        UA_Cold

    end

    properties(Dependent)
        powerNominalDischarge

        massRateCycCharge
        massRateCycDischarge
        massRateCycDischargeTurbine1
        massRateCycDischargeTurbine2

        capacityMassHot
        capacityMassCold
        massSaltHot1Ratio  % 热罐1质量比例
        massSaltHot2Ratio
        massSaltCold1Ratio
        massSaltCold2Ratio
        
        T2Charge
        T5Charge
        T6Charge
        T3Charge
        T2Discharge
        T5Discharge
        T6Discharge
        T3Discharge

        % 经济性参数
        investCostTotal
        investCostTurbine
        investCostCompressor
        investCostHeatExchanger
        investCostMotor
        investCostgenerator
        investCostTank
        investCostMaterial
        % investCostCooler 固态系统需要，LBPTES不需要

    end

    methods
        function obj = BPTES(powerNominalCharge, chargingHours, UA_Reg, UA_Hot, UA_Cold, T1Charge, T4Charge, T1Discharge, T4Discharge, pressureRatio, capacityMassGain)
            obj.powerNominalCharge = powerNominalCharge;
            obj.powerCharge = 1;
            obj.capacityEnergy = 3600 * chargingHours .* powerNominalCharge;
            obj.capacityMassGain = capacityMassGain;
            obj.UA_Reg = UA_Reg;
            obj.UA_Hot = UA_Hot;
            obj.UA_Cold = UA_Cold;
            obj.T1Charge = T1Charge;
            obj.T4Charge = T4Charge;
            obj.T1Discharge = T1Discharge;
            obj.T4Discharge = T4Discharge;
            obj.pressureRatio = pressureRatio;

            obj.massSaltHot1 =  0.1001.* obj.capacityEnergy ./ (obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot .* (1 + obj.capacityMassGain); % 热罐1质量
            obj.massSaltHot2 = 0.8999 .* obj.capacityEnergy ./ (obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot .* (1 + obj.capacityMassGain);
            obj.massSaltCold1 = 0.8999 .* obj.capacityEnergy ./ (obj.T1SaltCold - obj.T2SaltCold) ./ obj.c_pSaltCold .* (1 - obj.capacityMassGain);
            obj.massSaltCold2 = 0.1001 .* obj.capacityEnergy ./ (obj.T1SaltCold - obj.T2SaltCold) ./ obj.c_pSaltCold .* (1 - obj.capacityMassGain);
        end

        function value = get.powerNominalDischarge(obj)
            m_ch = obj.powerNominalCharge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            );

            m_dis = m_ch.*sqrt(obj.T2Charge./obj.T4Discharge);
            m_dis1 = m_ch.*sqrt(obj.T4Charge./obj.T1Discharge);
            m_dis2 = m_dis - m_dis1;
            % 计算单位质量工质吸热量与做功量
            heatIn = obj.c_pCycle .* (obj.T1Discharge - obj.T6Discharge); % 这里用的
            heatOut = obj.c_pCycle .* (obj.T3Discharge - obj.T4Discharge);
            workRate = heatIn - heatOut; % 单位质量工质做功功率

            value = ((m_dis1./m_dis).*obj.efficiencyGearBox + (m_dis2./m_dis)).*workRate.*m_dis.*obj.efficiencyEl;

        end

        function value = get.massRateCycCharge(obj)
            value = obj.powerCharge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            ); % 这里要求powerCharge初始化
        end

        function value = get.massRateCycDischarge(obj)
            value = obj.massRateCycCharge.*sqrt(obj.T2Charge./obj.T4Discharge); % > Charge，压缩机关联式
            % value = obj.massRateCycCharge.*sqrt(obj.T5Charge./obj.T1Discharge); % < Charge，透平关联式
        end

        function value = get.massRateCycDischargeTurbine1(obj)
            value = obj.massRateCycCharge .* sqrt(obj.T4Charge ./ obj.T1Discharge);
        end

        function value = get.massRateCycDischargeTurbine2(obj)
            value = obj.massRateCycDischarge - obj.massRateCycDischargeTurbine1;
        end

        function value = get.capacityMassHot(obj)
            value = obj.capacityEnergy ./ (obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot .* (1 + obj.capacityMassGain);
        end

        function value = get.capacityMassCold(obj)
            value = obj.capacityEnergy ./ (obj.T1SaltCold - obj.T2SaltCold) ./ obj.c_pSaltCold .* (1 - obj.capacityMassGain);
        end

        function value = get.massSaltHot1Ratio(obj)
            value = obj.massSaltHot1 ./ obj.capacityMassHot;
        end

        function value = get.massSaltHot2Ratio(obj)
            value = obj.massSaltHot2 ./ obj.capacityMassHot;
        end

        function value = get.massSaltCold1Ratio(obj)
            value = obj.massSaltCold1 ./ obj.capacityMassCold;
        end

        function value = get.massSaltCold2Ratio(obj)
            value = obj.massSaltCold2 ./ obj.capacityMassCold;
        end

        function value = get.T5Discharge(obj)
            T5_ = obj.T4Discharge .* (obj.pressureRatio .^ ((obj.gamma-1)./obj.gamma));
            value = obj.T4Discharge + (T5_ - obj.T4Discharge) ./ obj.efficiencyCompressor;
        end

        function value = get.T2Discharge(obj)
            T2_ = obj.T1Discharge ./ (obj.pressureRatio .^ ((obj.gamma-1)./obj.gamma));
            value = obj.T1Discharge - (obj.T1Discharge - T2_) .* obj.efficiencyTurbine;
        end

        function value = get.T2Charge(obj)
            value = obj.efficiencyCompressor ./ (obj.efficiencyCompressor + obj.pressureRatio.^ ((obj.gamma-1)./obj.gamma) -1) .* obj.T1Charge;
        end

        function value = get.T5Charge(obj)
            value = obj.T4Charge ./ (1-obj.efficiencyTurbine+obj.efficiencyTurbine./(obj.pressureRatio .^ ((obj.gamma-1)./obj.gamma)));
        end

        function value = get.T6Discharge(obj)
            value = (obj.T2Discharge + obj.massRateCycDischarge.*obj.c_pCycle./obj.UA_Reg./obj.f_Reg.*obj.T5Discharge) ./ (1 + obj.massRateCycDischarge.*obj.c_pCycle./obj.UA_Reg./obj.f_Reg);
            value = min(value, obj.T2Discharge - obj.pinchpointReg - obj.tol);
        end

        function value = get.T3Discharge(obj)
            value = obj.T2Discharge - (obj.T6Discharge - obj.T5Discharge);
            value = max(value, obj.T5Discharge + obj.pinchpointReg + obj.tol);
        end

        function value = get.T6Charge(obj)
            value = (obj.massRateCycCharge.*obj.c_pCycle./obj.UA_Reg./obj.f_Reg.* obj.T5Charge - obj.T2Charge) ./ (obj.massRateCycCharge.*obj.c_pCycle./obj.UA_Reg./obj.f_Reg - 1);
            value = max(value, obj.T2Charge + obj.pinchpointReg + obj.tol);
        end

        function value = get.T3Charge(obj)
            value = obj.T2Charge + obj.T5Charge - obj.T6Charge;
            value = min(value, obj.T5Charge - obj.pinchpointReg - obj.tol);
        end
        

        function powerReal = chargeORdischarge(obj, powerAim)
            if obj.massSaltHot1 < obj.capacityBound(1)*obj.capacityMassHot || ...
                obj.massSaltHot1 > obj.capacityBound(2)*obj.capacityMassHot || ...
                obj.massSaltHot2 < obj.capacityBound(1)*obj.capacityMassHot || ...
                obj.massSaltHot2 > obj.capacityBound(2)*obj.capacityMassHot || ...
                obj.massSaltCold1 < obj.capacityBound(1)*obj.capacityMassCold || ...
                obj.massSaltCold1 > obj.capacityBound(2)*obj.capacityMassCold || ...
                obj.massSaltCold2 < obj.capacityBound(1)*obj.capacityMassCold || ...
                obj.massSaltCold2 > obj.capacityBound(2)*obj.capacityMassCold

                disp('储热罐容量已经出现问题，无法进行充放电');
                disp(obj);
                % keyboard

            end

            if powerAim > obj.tol
                powerReal = obj.discharge(powerAim);
            elseif powerAim < - obj.tol
                powerReal = - obj.charge(powerAim);
            else
                powerReal = 0;
            end
        end

        function powerDischarge = discharge(obj, powerAim)
            powerAim = abs(powerAim);
            
            powerDischarge = min(powerAim, obj.powerNominalDischarge);

            % 计算单位质量工质吸热量与做功量
            heatIn = obj.c_pCycle .* (obj.T1Discharge - obj.T6Discharge); % 这里用的是回热温度，其实用无回热循环的温度/热量差也一样
            heatOut = obj.c_pCycle .* (obj.T3Discharge - obj.T4Discharge);
            workRate = heatIn - heatOut; % 单位质量工质做功功率

            for ii = 1:obj.maxIter+1
                powerDischarge = min(powerAim, obj.powerNominalDischarge) * (1 - (ii-1) / obj.maxIter);

                if ii == obj.maxIter
                    powerDischarge = obj.tol;
                end

                massRateCyc = powerDischarge ./ obj.efficiencyEl ./ workRate ./ ...
                ((obj.massRateCycDischargeTurbine1./obj.massRateCycDischarge).*obj.efficiencyGearBox + (obj.massRateCycDischargeTurbine2./obj.massRateCycDischarge));

                massHot1 = obj.massSaltHot1 - heatIn * massRateCyc ./(obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot * obj.timeStep;
                massHot2 = obj.massSaltHot2 + heatIn * massRateCyc ./(obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot * obj.timeStep;
                massCold1 = obj.massSaltCold1 + heatOut * massRateCyc ./(obj.T1SaltCold  - obj.T2SaltCold) ./ obj.c_pSaltCold * obj.timeStep;
                massCold2 = obj.massSaltCold2 - heatOut * massRateCyc ./(obj.T1SaltCold  - obj.T2SaltCold) ./ obj.c_pSaltCold * obj.timeStep;

                if massHot1 > obj.capacityBound(1)*obj.capacityMassHot && massHot1 < obj.capacityBound(2)*obj.capacityMassHot && ...
                    massHot2 > obj.capacityBound(1)*obj.capacityMassHot && massHot2 < obj.capacityBound(2)*obj.capacityMassHot && ...
                    massCold1 > obj.capacityBound(1)*obj.capacityMassCold && massCold1 < obj.capacityBound(2)*obj.capacityMassCold && ...
                    massCold2 > obj.capacityBound(1)*obj.capacityMassCold && massCold2 < obj.capacityBound(2)*obj.capacityMassCold
                    % massCold1 > obj.capacityBound(1)*obj.capacityMassCold  && ...
                    % massCold2 < obj.capacityBound(2)*obj.capacityMassCold % 考虑让冷罐可以选择不接入
                    break;
                end

                if ii == obj.maxIter+1
                    disp('放电过程达到最大迭代数，储热罐情况如下：');
                    disp(obj);
                    % keyboard;
                end
            end

            if ii < obj.maxIter+1
                obj.massSaltHot1 = massHot1;
                obj.massSaltHot2 = massHot2;
                obj.massSaltCold1 = massCold1;
                obj.massSaltCold2 = massCold2;
                % % 考虑冷罐可以不接入时
                % if massCold1 > obj.capacityBound(1)*obj.capacityMassCold && massCold1 < obj.capacityBound(2)*obj.capacityMassCold && ...
                %     massCold2 > obj.capacityBound(1)*obj.capacityMassCold && massCold2 < obj.capacityBound(2)*obj.capacityMassCold
                %     obj.massSaltCold1 = massCold1;
                %     obj.massSaltCold2 = massCold2; % 只有满足约束时，才改变冷罐容量
                % else
                %     obj.massSaltCold1 = obj.capacityBound(2)*obj.capacityMassCold - 1;
                %     obj.massSaltCold2 = obj.capacityBound(1)*obj.capacityMassCold + 1; % 放电时冷罐容量不变
                % end
                obj.powerDischarge = powerDischarge;
            end
        end

        function powerCharge = charge(obj, powerAim)
            powerAim = abs(powerAim);
            
            powerCharge = min(powerAim, obj.powerNominalCharge);

            % 计算单位质量工质吸热量与做功量
            heatIn = obj.c_pCycle .* (obj.T3Charge - obj.T4Charge); % 这里用的是回热温度，其实用无回热循环的温度/热量差也一样
            heatOut = obj.c_pCycle .* (obj.T1Charge - obj.T6Charge);
            workRate = heatOut - heatIn; % 单位质量工质做功功率

            for ii = 1:obj.maxIter+1
                powerCharge = min(powerAim, obj.powerNominalCharge) * (1 - (ii-1) / obj.maxIter);
                if ii == obj.maxIter
                    powerCharge = obj.tol;
                end

                massRateCyc = powerCharge ./ workRate .* obj.efficiencyEl .* obj.efficiencyGearBox;

                massHot1 = obj.massSaltHot1 + heatOut * massRateCyc ./(obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot * obj.timeStep;
                massHot2 = obj.massSaltHot2 - heatOut * massRateCyc ./(obj.T1SaltHot - obj.T2SaltHot) ./ obj.c_pSaltHot * obj.timeStep;
                massCold1 = obj.massSaltCold1 - heatIn * massRateCyc ./(obj.T1SaltCold - obj.T2SaltCold) ./ obj.c_pSaltCold * obj.timeStep;
                massCold2 = obj.massSaltCold2 + heatIn * massRateCyc ./(obj.T1SaltCold - obj.T2SaltCold) ./ obj.c_pSaltCold * obj.timeStep;

                if massHot1 > obj.capacityBound(1)*obj.capacityMassHot && massHot1 < obj.capacityBound(2)*obj.capacityMassHot && ...
                    massHot2 > obj.capacityBound(1)*obj.capacityMassHot && massHot2 < obj.capacityBound(2)*obj.capacityMassHot && ...
                    massCold1 > obj.capacityBound(1)*obj.capacityMassCold && massCold1 < obj.capacityBound(2)*obj.capacityMassCold && ...
                    massCold2 > obj.capacityBound(1)*obj.capacityMassCold && massCold2 < obj.capacityBound(2)*obj.capacityMassCold
                    % massCold1 < obj.capacityBound(2)*obj.capacityMassCold && ...
                    % massCold2 > obj.capacityBound(1)*obj.capacityMassCold
                    % 此处暂时考虑让冷罐超限
                    break;
                end

                if ii == obj.maxIter+1
                    disp('充电过程达到最大迭代数，储热罐情况如下：');
                    disp(obj);
                    % keyboard;
                end
            end

            if ii < obj.maxIter+1
                obj.massSaltHot1 = massHot1;
                obj.massSaltHot2 = massHot2;
                obj.massSaltCold1 = massCold1;
                obj.massSaltCold2 = massCold2;
                % % 考虑冷罐可以不接入时
                % if  massCold1 > obj.capacityBound(1)*obj.capacityMassCold && massCold1 < obj.capacityBound(2)*obj.capacityMassCold && ...
                %     massCold2 > obj.capacityBound(1)*obj.capacityMassCold && massCold2 < obj.capacityBound(2)*obj.capacityMassCold
                %     obj.massSaltCold1 = massCold1;
                %     obj.massSaltCold2 = massCold2; % 只有满足约束时，才改变冷罐容量
                %     else
                %     obj.massSaltCold1 = obj.capacityBound(1)*obj.capacityMassCold + 1;
                %     obj.massSaltCold2 = obj.capacityBound(2)*obj.capacityMassCold - 1; % 充电时冷罐容量不变
                % end
                obj.powerCharge = powerCharge;
            end

        end

        % 经济性
        function value = get.investCostTotal(obj)
            value = obj.investCostTurbine + obj.investCostCompressor + obj.investCostHeatExchanger + obj.investCostMotor + obj.investCostgenerator + obj.investCostTank + obj.investCostMaterial;
        end

        function value = get.investCostTurbine(obj)
            mass0 = obj.powerNominalDischarge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            ); % 单位kg/s
            mass0Charge = obj.powerNominalCharge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            ); % 单位kg/s
            mass1 = mass0.*sqrt(obj.T4Charge ./ obj.T1Discharge);
            cost1 = 1.051*266.3*max(mass1,mass0Charge).*log(obj.pressureRatio)./(0.94-obj.efficiencyTurbine);
            cost1 = cost1 .* (1 + exp(0.036*max(obj.T1Discharge, obj.T5Charge)-1.207*54.4));
            cost1 = cost1 .* (obj.pressureMin/obj.pressureAtm)^(-1) .*max(obj.T2Discharge, obj.T4Charge)./(300); % 理想气体假设
            mass2 = mass0 - mass1;
            cost2 = 1.051*266.3*mass2.*log(obj.pressureRatio)./(0.94-obj.efficiencyTurbine);
            cost2 = cost2 .* (1 + exp(0.036*obj.T1Discharge-1.207*54.4));
            cost2 = cost2 .* (obj.pressureMin/obj.pressureAtm)^(-1) .*max(obj.T2Discharge, obj.T4Charge)./(300); % 理想气体假设
            value = cost1 + cost2;
            CEPCI1995 = 381.1;
            CEPCI2022 = 813.0;
            value = value .* (CEPCI2022 / CEPCI1995);
        end

        function value = get.investCostCompressor(obj)
            fm2 = 2; % 与材料有关
            power0 = max(obj.powerNominalCharge, obj.powerNominalDischarge);
            mass0 = power0 ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            ); % 单位kg/s
            value = 1.051 * 39.5 * mass0 / (0.92 - obj.efficiencyCompressor) * obj.pressureRatio * log(obj.pressureRatio);
            % 假设入口压强等于大气压，用温度比例来计算密度比例
            value = value .* (obj.pressureMin/obj.pressureAtm)^(-1) .*max(obj.T4Discharge, obj.T2Charge)./(300); % 理想气体假设
            CEPCI1995 = 381.1;
            CEPCI2022 = 813.0;
            value = value .* (CEPCI2022 / CEPCI1995);
        end
 
        function value = get.investCostMotor(obj)
            % mass0Discharge = obj.powerNominalDischarge ./ ...
            %     (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Discharge-obj.T2Discharge) ...
            %     - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Discharge-obj.T4Discharge) ...
            % ); % 单位kg/s
            % mass0Charge = obj.powerNominalCharge ./ ...
            %     (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
            %     - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            % ); % 单位kg/s

            % power0 = max(obj.c_pCycle.*(obj.T5Discharge-obj.T4Discharge).*mass0Discharge./obj.efficiencyGearBox./obj.efficiencyEl, obj.c_pCycle.*(obj.T1Charge-obj.T2Charge).*mass0Charge./obj.efficiencyGearBox./obj.efficiencyEl);

            % value = 399400.*(power0./1e3).^0.61*1.4; % 1.4为齿轮箱修正系数
            % 以上为不对压缩机与透平用轴连接，直接利用发电机与电动机的算法
            % 用轴连接更为合理

            value = 399400.*(obj.powerNominalCharge./1e6).^0.61*1.4;

            CEPCI2017 = 567.5;
            CEPCI2022 = 813.0;            
            value = value .* (CEPCI2022 / CEPCI2017);
        end

        function value = get.investCostgenerator(obj)
            mass0 = obj.powerNominalDischarge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Discharge-obj.T2Discharge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Discharge-obj.T4Discharge) ...
            ); % 单位kg/s
            mass0Charge = obj.powerNominalCharge ./ ...
                (obj.c_pCycle./obj.efficiencyEl./obj.efficiencyGearBox.*(obj.T1Charge-obj.T2Charge) ...
                - obj.c_pCycle.*obj.efficiencyEl.*obj.efficiencyGearBox.*(obj.T5Charge-obj.T4Charge) ...
            ); % 单位kg/s
            mass1 = mass0Charge.*sqrt(obj.T4Charge ./ obj.T1Discharge);
            mass2 = mass0 - mass1;
            % mass1Charge = mass0Charge;
            % power1 = max(obj.c_pCycle.*(obj.T1Discharge-obj.T2Discharge).*obj.efficiencyGearBox.*obj.efficiencyEl.*mass1, obj.c_pCycle.*(obj.T5Charge-obj.T4Charge).*obj.efficiencyGearBox.*obj.efficiencyEl.*mass1Charge);

            % power2 = mass2.*obj.c_pCycle.*(obj.T1Discharge-obj.T2Discharge).*obj.efficiencyEl;

            % cost1 = 108900.*(power1./1e3).^0.55.*1.4; % 1.4为齿轮箱修正系数
            % cost2 = 108900.*(power2./1e3).^0.55;
            % 以上为不对压缩机与透平用轴连接，直接利用发电机与电动机的算法
            % 用轴连接更为合理

            cost1 = 108900.*(mass1./mass0.*obj.powerNominalDischarge./1e6).^0.55.*1.4; % 1.4为齿轮箱修正系数
            cost2 = 108900.*(mass2./mass0.*obj.powerNominalDischarge./1e6).^0.55;

            value = cost1 + cost2;
            CEPCI2017 = 567.5;
            CEPCI2022 = 813.0;            
            value = value .* (CEPCI2022 / CEPCI2017);
        end

        function value = get.investCostHeatExchanger(obj)
            value = 49.45.*(obj.UA_Cold).^0.75 + 49.45.*(obj.UA_Hot).^0.75 + 49.45.*(obj.UA_Reg).^0.75;
            CEPCI2017 = 567.5;
            CEPCI2022 = 813.0;            
            value = value .* (CEPCI2022 / CEPCI2017);
        end

        function value = densitySalt(obj,T)
            if T > 273.15 + 150
                value = 2090 - 0.636 * (T - 273.15);         % 熔融盐密度，kg/m^3
            elseif T < 273.15+50
                value = 918.1093 - 0.8807 * (T - 273.15);
            else
                error('储热介质温度出现在奇怪的区间');
            end
        end

        function value = get.investCostTank(obj)
            f_m = 1;
            V_tank = [obj.capacityMassHot./obj.densitySalt(obj.T1SaltHot), obj.capacityMassHot./obj.densitySalt(obj.T2SaltHot), obj.capacityMassCold./obj.densitySalt(obj.T1SaltCold), obj.capacityMassCold./obj.densitySalt(obj.T2SaltCold)];
            value = sum(f_m.* (170.5 .* V_tank + 59560));
            CEPCI1990 = 357.6;
            CEPCI2022 = 813.0;            
            value = value .* (CEPCI2022 ./ CEPCI1990);
        end

        function value = get.investCostMaterial(obj)
            value = 0.61 * obj.capacityMassCold + 0.54 * obj.capacityMassHot;
            CEPCI2016 = 541.7;
            CEPCI2022 = 813.0;            
            value = value .* (CEPCI2022 / CEPCI2016);
        end







    end
end
